In [ ]:
## Alyssa Felix-Arreola

Assignment 4: Street Networks & Web Scraping¶

Part 1: Visualizing crash data in Philadelphia

In this section, you will use osmnx to analyze the crash incidence in Center City.

Part 2: Scraping Craigslist

In this section, you will use Selenium and BeautifulSoup to scrape data for hundreds of apartments from Philadelphia's Craigslist portal.

In [132]:
import pandas as pd
import geopandas as gpd
import numpy as np
import holoviews as hv
import panel as pn
from matplotlib import pyplot as plt
In [133]:
# Show all columns in dataframes
pd.options.display.max_columns = 999
In [134]:
# Hide warnings due to issue in shapely package 
# See: https://github.com/shapely/shapely/issues/1345
np.seterr(invalid="ignore");

Part 1: Visualizing crash data in Philadelphia¶

1.1 Load the geometry for the region being analyzed¶

We'll analyze crashes in the "Central" planning district in Philadelphia, a rough approximation for Center City. Planning districts can be loaded from Open Data Philly. Read the data into a GeoDataFrame using the following link:

http://data.phl.opendata.arcgis.com/datasets/0960ea0f38f44146bb562f2b212075aa_0.geojson

Select the "Central" district and extract the geometry polygon for only this district. After this part, you should have a polygon variable of type shapely.geometry.polygon.Polygon.

In [135]:
planningdistricts = gpd.read_file('~/miniforge3/envs/musa-550-fall-2023/assignment-4-alyssafelixa-main/Planning_Districts.geojson')
central= planningdistricts.query('DIST_NAME == "Central"')
central = central.squeeze().geometry

1.2 Get the street network graph¶

Use OSMnx to create a network graph (of type 'drive') from your polygon boundary in 1.1.

In [136]:
import osmnx as ox
In [137]:
cstreets = ox.graph_from_polygon(central, network_type='drive')

ox.plot_graph(cstreets, node_size=0, edge_color='b', figsize=(10, 10))
No description has been provided for this image
Out[137]:
(<Figure size 1000x1000 with 1 Axes>, <Axes: >)

1.3 Convert your network graph edges to a GeoDataFrame¶

Use OSMnx to create a GeoDataFrame of the network edges in the graph object from part 1.2. The GeoDataFrame should contain the edges but not the nodes from the network.

In [138]:
cc_edges = ox.graph_to_gdfs(cstreets, edges=True, nodes=False)
cc_edges.head()
Out[138]:
osmid oneway lanes name highway reversed length geometry maxspeed bridge ref tunnel width service access junction
u v key
109727439 109911666 0 132508434 True 1 Bainbridge Street residential False 44.347 LINESTRING (-75.17104 39.94345, -75.17060 39.9... NaN NaN NaN NaN NaN NaN NaN NaN
109727448 109727439 0 12109011 True NaN South Colorado Street residential False 109.496 LINESTRING (-75.17125 39.94248, -75.17120 39.9... NaN NaN NaN NaN NaN NaN NaN NaN
110034229 0 12159387 True NaN Fitzwater Street residential False 91.353 LINESTRING (-75.17125 39.94248, -75.17137 39.9... NaN NaN NaN NaN NaN NaN NaN NaN
109727507 110024052 0 193364514 True NaN Carpenter Street residential False 53.208 LINESTRING (-75.17196 39.93973, -75.17134 39.9... NaN NaN NaN NaN NaN NaN NaN NaN
109728761 110274344 0 672312336 True NaN Brown Street residential False 58.270 LINESTRING (-75.17317 39.96951, -75.17250 39.9... 25 mph NaN NaN NaN NaN NaN NaN NaN
In [139]:
ax = cc_edges.to_crs(epsg=3857).plot(color="gray")
boundary = gpd.GeoSeries([central], crs="EPSG:4326")
boundary.to_crs(epsg=3857).plot(
    ax=ax, facecolor="none", edgecolor="red", linewidth=3, zorder=2
)

ax.set_axis_off()
No description has been provided for this image

1.4 Load PennDOT crash data¶

Data for crashes (of all types) for 2020, 2021, and 2022 in Philadelphia County is available at the following path:

./data/CRASH_PHILADELPHIA_XXXX.csv

You should see three separate files in the data/ folder. Use pandas to read each of the CSV files, and combine them into a single dataframe using pd.concat().

The data was downloaded for Philadelphia County from here.

In [140]:
crash2020 = pd.read_csv('~/miniforge3/envs/musa-550-fall-2023/assignment-4-alyssafelixa-main/data/CRASH_PHILADELPHIA_2020.csv')
crash2021 = pd.read_csv('~/miniforge3/envs/musa-550-fall-2023/assignment-4-alyssafelixa-main/data/CRASH_PHILADELPHIA_2021.csv')
crash2022 = pd.read_csv('~/miniforge3/envs/musa-550-fall-2023/assignment-4-alyssafelixa-main/data/CRASH_PHILADELPHIA_2022.csv')
crash2020.head()
Out[140]:
CRN ARRIVAL_TM AUTOMOBILE_COUNT BELTED_DEATH_COUNT BELTED_SUSP_SERIOUS_INJ_COUNT BICYCLE_COUNT BICYCLE_DEATH_COUNT BICYCLE_SUSP_SERIOUS_INJ_COUNT BUS_COUNT CHLDPAS_DEATH_COUNT CHLDPAS_SUSP_SERIOUS_INJ_COUNT COLLISION_TYPE COMM_VEH_COUNT CONS_ZONE_SPD_LIM COUNTY CRASH_MONTH CRASH_YEAR DAY_OF_WEEK DEC_LAT DEC_LONG DISPATCH_TM DISTRICT DRIVER_COUNT_16YR DRIVER_COUNT_17YR DRIVER_COUNT_18YR DRIVER_COUNT_19YR DRIVER_COUNT_20YR DRIVER_COUNT_50_64YR DRIVER_COUNT_65_74YR DRIVER_COUNT_75PLUS EST_HRS_CLOSED FATAL_COUNT HEAVY_TRUCK_COUNT HORSE_BUGGY_COUNT HOUR_OF_DAY ILLUMINATION INJURY_COUNT INTERSECT_TYPE INTERSECTION_RELATED LANE_CLOSED LATITUDE LN_CLOSE_DIR LOCATION_TYPE LONGITUDE MAX_SEVERITY_LEVEL MCYCLE_DEATH_COUNT MCYCLE_SUSP_SERIOUS_INJ_COUNT MOTORCYCLE_COUNT MUNICIPALITY NONMOTR_COUNT NONMOTR_DEATH_COUNT NONMOTR_SUSP_SERIOUS_INJ_COUNT NTFY_HIWY_MAINT PED_COUNT PED_DEATH_COUNT PED_SUSP_SERIOUS_INJ_COUNT PERSON_COUNT POLICE_AGCY POSSIBLE_INJ_COUNT RDWY_SURF_TYPE_CD RELATION_TO_ROAD ROAD_CONDITION ROADWAY_CLEARED SCH_BUS_IND SCH_ZONE_IND SECONDARY_CRASH SMALL_TRUCK_COUNT SPEC_JURIS_CD SUSP_MINOR_INJ_COUNT SUSP_SERIOUS_INJ_COUNT SUV_COUNT TCD_FUNC_CD TCD_TYPE TFC_DETOUR_IND TIME_OF_DAY TOT_INJ_COUNT TOTAL_UNITS UNB_DEATH_COUNT UNB_SUSP_SERIOUS_INJ_COUNT UNBELTED_OCC_COUNT UNK_INJ_DEG_COUNT UNK_INJ_PER_COUNT URBAN_AREA URBAN_RURAL VAN_COUNT VEHICLE_COUNT WEATHER1 WEATHER2 WORK_ZONE_IND WORK_ZONE_LOC WORK_ZONE_TYPE WORKERS_PRES WZ_CLOSE_DETOUR WZ_FLAGGER WZ_LAW_OFFCR_IND WZ_LN_CLOSURE WZ_MOVING WZ_OTHER WZ_SHLDER_MDN WZ_WORKERS_INJ_KILLED
0 2020036588 1349.0 1 0 0 0 0 0 0 0 0 1 0 NaN 67 3 2020 2 39.9601 -75.1794 1343.0 6 0 0 0 0 0 0 0 0 NaN 0 0 0.0 13 1 0 0 NaN 1 39 57:36.245 4.0 0 75 10:45.819 0 0 0 0 67301 0 0 0 N 0 0 0 3 68K01 0 NaN 1 9 NaN N N NaN 1 NaN 0 0 0 0 0 N 1332 0 2 0 0 0 0 0 3 4 0 2 7 NaN N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
1 2020036617 1842.0 1 0 0 0 0 0 0 0 0 7 0 NaN 67 4 2020 1 39.9809 -75.2065 1840.0 6 0 0 0 0 0 0 0 0 NaN 0 0 0.0 18 1 0 0 NaN 1 39 58:51.367 3.0 0 75 12:23.366 0 0 0 0 67301 0 0 0 N 0 0 0 1 68K01 0 NaN 2 9 NaN N N NaN 0 NaN 0 0 0 0 0 N 1838 0 1 0 0 0 0 0 3 4 0 1 7 7.0 N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2 2020035717 2000.0 1 0 0 0 0 0 0 0 0 2 0 NaN 67 4 2020 1 39.9269 -75.1691 2000.0 6 0 0 0 0 0 1 0 0 NaN 0 0 0.0 14 1 1 1 NaN 9 39 55:36.660 NaN 0 75 10:08.677 3 0 0 0 67301 0 0 0 N 0 0 0 2 67301 0 NaN 1 9 NaN N N NaN 0 NaN 1 0 1 3 3 U 1457 1 2 0 0 0 0 0 3 4 0 2 7 4.0 N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
3 2020034378 1139.0 2 0 0 0 0 0 0 0 0 1 0 NaN 67 4 2020 4 39.9237 -75.1924 1131.0 6 0 0 0 0 0 2 0 0 NaN 0 0 0.0 11 1 1 0 NaN 0 39 55:25.183 NaN 0 75 11:32.758 3 0 0 0 67301 0 0 0 N 0 0 0 3 68K01 0 NaN 1 1 NaN N N NaN 1 NaN 1 0 0 0 0 NaN 1128 1 3 0 0 0 0 0 3 4 0 3 3 NaN N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
4 2020025511 345.0 1 0 0 0 0 0 0 0 0 7 0 NaN 67 3 2020 1 39.8826 -75.2450 329.0 6 0 0 0 0 0 0 0 0 NaN 0 0 0.0 3 3 2 0 NaN 0 39 52:57.717 NaN 0 75 14:41.931 3 0 0 0 67301 0 0 0 Y 0 0 0 2 68K01 0 NaN 4 1 NaN N N NaN 0 NaN 2 0 0 0 0 NaN 328 2 1 0 0 0 0 0 3 4 0 1 3 3.0 N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
In [141]:
crashes = pd.concat([crash2020, crash2021, crash2022], ignore_index=True)
crashes.head()
Out[141]:
CRN ARRIVAL_TM AUTOMOBILE_COUNT BELTED_DEATH_COUNT BELTED_SUSP_SERIOUS_INJ_COUNT BICYCLE_COUNT BICYCLE_DEATH_COUNT BICYCLE_SUSP_SERIOUS_INJ_COUNT BUS_COUNT CHLDPAS_DEATH_COUNT CHLDPAS_SUSP_SERIOUS_INJ_COUNT COLLISION_TYPE COMM_VEH_COUNT CONS_ZONE_SPD_LIM COUNTY CRASH_MONTH CRASH_YEAR DAY_OF_WEEK DEC_LAT DEC_LONG DISPATCH_TM DISTRICT DRIVER_COUNT_16YR DRIVER_COUNT_17YR DRIVER_COUNT_18YR DRIVER_COUNT_19YR DRIVER_COUNT_20YR DRIVER_COUNT_50_64YR DRIVER_COUNT_65_74YR DRIVER_COUNT_75PLUS EST_HRS_CLOSED FATAL_COUNT HEAVY_TRUCK_COUNT HORSE_BUGGY_COUNT HOUR_OF_DAY ILLUMINATION INJURY_COUNT INTERSECT_TYPE INTERSECTION_RELATED LANE_CLOSED LATITUDE LN_CLOSE_DIR LOCATION_TYPE LONGITUDE MAX_SEVERITY_LEVEL MCYCLE_DEATH_COUNT MCYCLE_SUSP_SERIOUS_INJ_COUNT MOTORCYCLE_COUNT MUNICIPALITY NONMOTR_COUNT NONMOTR_DEATH_COUNT NONMOTR_SUSP_SERIOUS_INJ_COUNT NTFY_HIWY_MAINT PED_COUNT PED_DEATH_COUNT PED_SUSP_SERIOUS_INJ_COUNT PERSON_COUNT POLICE_AGCY POSSIBLE_INJ_COUNT RDWY_SURF_TYPE_CD RELATION_TO_ROAD ROAD_CONDITION ROADWAY_CLEARED SCH_BUS_IND SCH_ZONE_IND SECONDARY_CRASH SMALL_TRUCK_COUNT SPEC_JURIS_CD SUSP_MINOR_INJ_COUNT SUSP_SERIOUS_INJ_COUNT SUV_COUNT TCD_FUNC_CD TCD_TYPE TFC_DETOUR_IND TIME_OF_DAY TOT_INJ_COUNT TOTAL_UNITS UNB_DEATH_COUNT UNB_SUSP_SERIOUS_INJ_COUNT UNBELTED_OCC_COUNT UNK_INJ_DEG_COUNT UNK_INJ_PER_COUNT URBAN_AREA URBAN_RURAL VAN_COUNT VEHICLE_COUNT WEATHER1 WEATHER2 WORK_ZONE_IND WORK_ZONE_LOC WORK_ZONE_TYPE WORKERS_PRES WZ_CLOSE_DETOUR WZ_FLAGGER WZ_LAW_OFFCR_IND WZ_LN_CLOSURE WZ_MOVING WZ_OTHER WZ_SHLDER_MDN WZ_WORKERS_INJ_KILLED
0 2020036588 1349.0 1 0 0 0 0 0 0 0 0 1 0 NaN 67 3 2020 2 39.9601 -75.1794 1343.0 6 0 0 0 0 0 0 0 0 NaN 0 0 0.0 13 1 0 0 NaN 1 39 57:36.245 4.0 0 75 10:45.819 0 0 0 0 67301 0 0 0 N 0 0 0 3 68K01 0 NaN 1 9 NaN N N NaN 1 NaN 0 0 0 0 0 N 1332 0 2 0 0 0 0 0 3 4 0 2 7 NaN N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
1 2020036617 1842.0 1 0 0 0 0 0 0 0 0 7 0 NaN 67 4 2020 1 39.9809 -75.2065 1840.0 6 0 0 0 0 0 0 0 0 NaN 0 0 0.0 18 1 0 0 NaN 1 39 58:51.367 3.0 0 75 12:23.366 0 0 0 0 67301 0 0 0 N 0 0 0 1 68K01 0 NaN 2 9 NaN N N NaN 0 NaN 0 0 0 0 0 N 1838 0 1 0 0 0 0 0 3 4 0 1 7 7.0 N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2 2020035717 2000.0 1 0 0 0 0 0 0 0 0 2 0 NaN 67 4 2020 1 39.9269 -75.1691 2000.0 6 0 0 0 0 0 1 0 0 NaN 0 0 0.0 14 1 1 1 NaN 9 39 55:36.660 NaN 0 75 10:08.677 3 0 0 0 67301 0 0 0 N 0 0 0 2 67301 0 NaN 1 9 NaN N N NaN 0 NaN 1 0 1 3 3 U 1457 1 2 0 0 0 0 0 3 4 0 2 7 4.0 N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
3 2020034378 1139.0 2 0 0 0 0 0 0 0 0 1 0 NaN 67 4 2020 4 39.9237 -75.1924 1131.0 6 0 0 0 0 0 2 0 0 NaN 0 0 0.0 11 1 1 0 NaN 0 39 55:25.183 NaN 0 75 11:32.758 3 0 0 0 67301 0 0 0 N 0 0 0 3 68K01 0 NaN 1 1 NaN N N NaN 1 NaN 1 0 0 0 0 NaN 1128 1 3 0 0 0 0 0 3 4 0 3 3 NaN N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
4 2020025511 345.0 1 0 0 0 0 0 0 0 0 7 0 NaN 67 3 2020 1 39.8826 -75.2450 329.0 6 0 0 0 0 0 0 0 0 NaN 0 0 0.0 3 3 2 0 NaN 0 39 52:57.717 NaN 0 75 14:41.931 3 0 0 0 67301 0 0 0 Y 0 0 0 2 68K01 0 NaN 4 1 NaN N N NaN 0 NaN 2 0 0 0 0 NaN 328 2 1 0 0 0 0 0 3 4 0 1 3 3.0 N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

1.5 Convert the crash data to a GeoDataFrame¶

You will need to use the DEC_LAT and DEC_LONG columns for latitude and longitude.

The full data dictionary for the data is available here

In [142]:
geometry = gpd.points_from_xy(crashes['DEC_LONG'], crashes['DEC_LAT'])
crashesgdf = gpd.GeoDataFrame(crashes, geometry=geometry, crs="EPSG:4326")
crashesgdf.head()
Out[142]:
CRN ARRIVAL_TM AUTOMOBILE_COUNT BELTED_DEATH_COUNT BELTED_SUSP_SERIOUS_INJ_COUNT BICYCLE_COUNT BICYCLE_DEATH_COUNT BICYCLE_SUSP_SERIOUS_INJ_COUNT BUS_COUNT CHLDPAS_DEATH_COUNT CHLDPAS_SUSP_SERIOUS_INJ_COUNT COLLISION_TYPE COMM_VEH_COUNT CONS_ZONE_SPD_LIM COUNTY CRASH_MONTH CRASH_YEAR DAY_OF_WEEK DEC_LAT DEC_LONG DISPATCH_TM DISTRICT DRIVER_COUNT_16YR DRIVER_COUNT_17YR DRIVER_COUNT_18YR DRIVER_COUNT_19YR DRIVER_COUNT_20YR DRIVER_COUNT_50_64YR DRIVER_COUNT_65_74YR DRIVER_COUNT_75PLUS EST_HRS_CLOSED FATAL_COUNT HEAVY_TRUCK_COUNT HORSE_BUGGY_COUNT HOUR_OF_DAY ILLUMINATION INJURY_COUNT INTERSECT_TYPE INTERSECTION_RELATED LANE_CLOSED LATITUDE LN_CLOSE_DIR LOCATION_TYPE LONGITUDE MAX_SEVERITY_LEVEL MCYCLE_DEATH_COUNT MCYCLE_SUSP_SERIOUS_INJ_COUNT MOTORCYCLE_COUNT MUNICIPALITY NONMOTR_COUNT NONMOTR_DEATH_COUNT NONMOTR_SUSP_SERIOUS_INJ_COUNT NTFY_HIWY_MAINT PED_COUNT PED_DEATH_COUNT PED_SUSP_SERIOUS_INJ_COUNT PERSON_COUNT POLICE_AGCY POSSIBLE_INJ_COUNT RDWY_SURF_TYPE_CD RELATION_TO_ROAD ROAD_CONDITION ROADWAY_CLEARED SCH_BUS_IND SCH_ZONE_IND SECONDARY_CRASH SMALL_TRUCK_COUNT SPEC_JURIS_CD SUSP_MINOR_INJ_COUNT SUSP_SERIOUS_INJ_COUNT SUV_COUNT TCD_FUNC_CD TCD_TYPE TFC_DETOUR_IND TIME_OF_DAY TOT_INJ_COUNT TOTAL_UNITS UNB_DEATH_COUNT UNB_SUSP_SERIOUS_INJ_COUNT UNBELTED_OCC_COUNT UNK_INJ_DEG_COUNT UNK_INJ_PER_COUNT URBAN_AREA URBAN_RURAL VAN_COUNT VEHICLE_COUNT WEATHER1 WEATHER2 WORK_ZONE_IND WORK_ZONE_LOC WORK_ZONE_TYPE WORKERS_PRES WZ_CLOSE_DETOUR WZ_FLAGGER WZ_LAW_OFFCR_IND WZ_LN_CLOSURE WZ_MOVING WZ_OTHER WZ_SHLDER_MDN WZ_WORKERS_INJ_KILLED geometry
0 2020036588 1349.0 1 0 0 0 0 0 0 0 0 1 0 NaN 67 3 2020 2 39.9601 -75.1794 1343.0 6 0 0 0 0 0 0 0 0 NaN 0 0 0.0 13 1 0 0 NaN 1 39 57:36.245 4.0 0 75 10:45.819 0 0 0 0 67301 0 0 0 N 0 0 0 3 68K01 0 NaN 1 9 NaN N N NaN 1 NaN 0 0 0 0 0 N 1332 0 2 0 0 0 0 0 3 4 0 2 7 NaN N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN POINT (-75.17940 39.96010)
1 2020036617 1842.0 1 0 0 0 0 0 0 0 0 7 0 NaN 67 4 2020 1 39.9809 -75.2065 1840.0 6 0 0 0 0 0 0 0 0 NaN 0 0 0.0 18 1 0 0 NaN 1 39 58:51.367 3.0 0 75 12:23.366 0 0 0 0 67301 0 0 0 N 0 0 0 1 68K01 0 NaN 2 9 NaN N N NaN 0 NaN 0 0 0 0 0 N 1838 0 1 0 0 0 0 0 3 4 0 1 7 7.0 N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN POINT (-75.20650 39.98090)
2 2020035717 2000.0 1 0 0 0 0 0 0 0 0 2 0 NaN 67 4 2020 1 39.9269 -75.1691 2000.0 6 0 0 0 0 0 1 0 0 NaN 0 0 0.0 14 1 1 1 NaN 9 39 55:36.660 NaN 0 75 10:08.677 3 0 0 0 67301 0 0 0 N 0 0 0 2 67301 0 NaN 1 9 NaN N N NaN 0 NaN 1 0 1 3 3 U 1457 1 2 0 0 0 0 0 3 4 0 2 7 4.0 N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN POINT (-75.16910 39.92690)
3 2020034378 1139.0 2 0 0 0 0 0 0 0 0 1 0 NaN 67 4 2020 4 39.9237 -75.1924 1131.0 6 0 0 0 0 0 2 0 0 NaN 0 0 0.0 11 1 1 0 NaN 0 39 55:25.183 NaN 0 75 11:32.758 3 0 0 0 67301 0 0 0 N 0 0 0 3 68K01 0 NaN 1 1 NaN N N NaN 1 NaN 1 0 0 0 0 NaN 1128 1 3 0 0 0 0 0 3 4 0 3 3 NaN N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN POINT (-75.19240 39.92370)
4 2020025511 345.0 1 0 0 0 0 0 0 0 0 7 0 NaN 67 3 2020 1 39.8826 -75.2450 329.0 6 0 0 0 0 0 0 0 0 NaN 0 0 0.0 3 3 2 0 NaN 0 39 52:57.717 NaN 0 75 14:41.931 3 0 0 0 67301 0 0 0 Y 0 0 0 2 68K01 0 NaN 4 1 NaN N N NaN 0 NaN 2 0 0 0 0 NaN 328 2 1 0 0 0 0 0 3 4 0 1 3 3.0 N NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN POINT (-75.24500 39.88260)

1.6 Trim the crash data to Center City¶

  1. Get the boundary of the edges data frame (from part 1.3). Accessing the .geometry.unary_union.convex_hull property will give you a nice outer boundary region.
  2. Trim the crashes using the within() function of the crash GeoDataFrame to find which crashes are within the boundary.

There should be about 3,750 crashes within the Central district.

In [143]:
cc_edges = cc_edges.to_crs(crashesgdf.crs)
In [144]:
conhull = cc_edges.geometry.unary_union.convex_hull
crashes_boundary = crashesgdf[crashesgdf.within(conhull)]

len(crashes_boundary)
Out[144]:
3751

1.7 Re-project our data into an approriate CRS

We'll need to find the nearest edge (street) in our graph for each crash. To do this, osmnx will calculate the distance from each crash to the graph edges. For this calculation to be accurate, we need to convert from latitude/longitude

We'll convert the local state plane CRS for Philadelphia, EPSG=2272

Two steps:¶

  1. Project the graph object (G) using the ox.project_graph. Run ox.project_graph? to see the documentation for how to convert to a specific CRS.
  2. Project the crash data using the .to_crs() function.
In [145]:
cstreets = ox.project_graph(cstreets, to_crs='epsg:2272')
In [146]:
crashes_boundary = crashes_boundary.to_crs(epsg=2272)
cc_edges = cc_edges.to_crs(epsg=2272)

1.8 Find the nearest edge for each crash¶

See: ox.distance.nearest_edges(). It takes three arguments:

  • the network graph
  • the longitude of your crash data (the x attribute of the geometry column)
  • the latitude of your crash data (the y attribute of the geometry column)

You will get a numpy array with 3 columns that represent (u, v, key) where each u and v are the node IDs that the edge links together. We will ignore the key value for our analysis.

In [147]:
nearest_crash = ox.distance.nearest_edges(cstreets, crashes_boundary.geometry.x, crashes_boundary.geometry.y)

1.9 Calculate the total number of crashes per street¶

  1. Make a DataFrame from your data from part 1.7 with three columns, u, v, and key (we will only use the u and v columns)
  2. Group by u and v and calculate the size
  3. Reset the index and name your size() column as crash_count

After this step you should have a DataFrame with three columns: u, v, and crash_count.

In [148]:
nearestcrash = pd.DataFrame(nearest_crash, columns=['u', 'v', 'key'])
nearestcrash.head()
Out[148]:
u v key
0 8482829382 7065714513 0
1 110369693 109848091 0
2 110161185 2043885885 0
3 775424610 775424581 0
4 110416511 110417392 0
In [149]:
street_crashes = nearestcrash.groupby(['u', 'v', 'key']).size().reset_index(name='crash_count')

street_crashes.head()
Out[149]:
u v key crash_count
0 109729330 110216446 0 1
1 109729474 3425014859 0 2
2 109729486 110342146 0 4
3 109729673 11345699563 0 3
4 109729699 11345699562 0 3

1.10 Merge your edges GeoDataFrame and crash count DataFrame¶

You can use pandas to merge them on the u and v columns. This will associate the total crash count with each edge in the street network.

Tips:

  • Use a left merge where the first argument of the merge is the edges GeoDataFrame. This ensures no edges are removed during the merge.
  • Use the fillna(0) function to fill in missing crash count values with zero.
In [150]:
merged_edges = cc_edges.merge(street_crashes, how='left', on=['u', 'v'])
merged_edges['crash_count'] = merged_edges['crash_count'].fillna(0)
merged_edges.head()
Out[150]:
u v osmid oneway lanes name highway reversed length geometry maxspeed bridge ref tunnel width service access junction key crash_count
0 109727439 109911666 132508434 True 1 Bainbridge Street residential False 44.347 LINESTRING (2691520.089 232816.131, 2691645.19... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
1 109727448 109727439 12109011 True NaN South Colorado Street residential False 109.496 LINESTRING (2691472.521 232460.544, 2691484.46... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
2 109727448 110034229 12159387 True NaN Fitzwater Street residential False 91.353 LINESTRING (2691472.521 232460.544, 2691438.30... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
3 109727507 110024052 193364514 True NaN Carpenter Street residential False 53.208 LINESTRING (2691302.734 231454.473, 2691476.43... NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.0
4 109728761 110274344 672312336 True NaN Brown Street residential False 58.270 LINESTRING (2690646.511 242288.027, 2690836.44... 25 mph NaN NaN NaN NaN NaN NaN NaN NaN 0.0

1.11 Calculate a "Crash Index"¶

Let's calculate a "crash index" that provides a normalized measure of the crash frequency per street. To do this, we'll need to:

  1. Calculate the total crash count divided by the street length, using the length column
  2. Perform a log transformation of the crash/length variable — use numpy's log10() function
  3. Normalize the index from 0 to 1 (see the lecture notes for an example of this transformation)

Note: since the crash index involves a log transformation, you should only calculate the index for streets where the crash count is greater than zero.

After this step, you should have a new column in the data frame from 1.9 that includes a column called part 1.9.

In [154]:
edges_crash = merged_edges[merged_edges['crash_count'] > 0]
merged_edges.loc[edges_crash.index, 'crash_index'] = np.log10(valid_rows['crash_count'] / valid_rows['length'])
merged_edges['normalized_crash_index'] = (merged_edges['crash_index'] - merged_edges['crash_index'].min()) / (merged_edges['crash_index'].max() - merged_edges['crash_index'].min())
print(merged_edges.head())

           u          v      osmid  oneway lanes                   name  \
0  109727439  109911666  132508434    True     1      Bainbridge Street   
1  109727448  109727439   12109011    True   NaN  South Colorado Street   
2  109727448  110034229   12159387    True   NaN       Fitzwater Street   
3  109727507  110024052  193364514    True   NaN       Carpenter Street   
4  109728761  110274344  672312336    True   NaN           Brown Street   

       highway reversed   length  \
0  residential    False   44.347   
1  residential    False  109.496   
2  residential    False   91.353   
3  residential    False   53.208   
4  residential    False   58.270   

                                            geometry maxspeed bridge  ref  \
0  LINESTRING (2691520.089 232816.131, 2691645.19...      NaN    NaN  NaN   
1  LINESTRING (2691472.521 232460.544, 2691484.46...      NaN    NaN  NaN   
2  LINESTRING (2691472.521 232460.544, 2691438.30...      NaN    NaN  NaN   
3  LINESTRING (2691302.734 231454.473, 2691476.43...      NaN    NaN  NaN   
4  LINESTRING (2690646.511 242288.027, 2690836.44...   25 mph    NaN  NaN   

  tunnel width service access junction  key  crash_count  \
0    NaN   NaN     NaN    NaN      NaN  NaN          0.0   
1    NaN   NaN     NaN    NaN      NaN  NaN          0.0   
2    NaN   NaN     NaN    NaN      NaN  NaN          0.0   
3    NaN   NaN     NaN    NaN      NaN  NaN          0.0   
4    NaN   NaN     NaN    NaN      NaN  NaN          0.0   

   normalized_crash_index  crash_index  
0                     NaN          NaN  
1                     NaN          NaN  
2                     NaN          NaN  
3                     NaN          NaN  
4                     NaN          NaN

1.12 Plot a histogram of the crash index values¶

Use matplotlib's hist() function to plot the crash index values from the previous step.

You should see that the index values are Gaussian-distributed, providing justification for why we log-transformed!

In [155]:
import matplotlib.pyplot as plt

crash_index_values = merged_edges['normalized_crash_index'].dropna()

plt.hist(crash_index_values, bins=30, color='blue', edgecolor='black')
plt.title('Histogram of Normalized Crash Index')
plt.xlabel('Normalized Crash Index')
plt.ylabel('Frequency')
plt.show()
No description has been provided for this image

1.13 Plot an interactive map of the street networks, colored by the crash index¶

You can use GeoPandas to make an interactive Folium map, coloring the streets by the crash index column.

Tip: if you use the viridis color map, try using a "dark" tile set for better constrast of the colors.

In [156]:
merged_edges.explore(
    column="normalized_crash_index", 
    cmap="viridis", 
    tiles="CartoDB positron", 
)
Out[156]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Part 2: Scraping Craigslist¶

In this part, we'll be extracting information on apartments from Craigslist search results. You'll be using Selenium and BeautifulSoup to extract the relevant information from the HTML text.

For reference on CSS selectors, please see the notes from Week 6.

Primer: the Craigslist website URL¶

We'll start with the Philadelphia region. First we need to figure out how to submit a query to Craigslist. As with many websites, one way you can do this is simply by constructing the proper URL and sending it to Craigslist.

https://philadelphia.craigslist.org/search/apa?min_price=1&min_bedrooms=1&minSqft=1#search=1~gallery~0~0

There are three components to this URL.

  1. The base URL: http://philadelphia.craigslist.org/search/apa

  2. The user's search parameters: ?min_price=1&min_bedrooms=1&minSqft=1

We will send nonzero defaults for some parameters (bedrooms, size, price) in order to exclude results that have empty values for these parameters.

  1. The URL hash: #search=1~gallery~0~0

As we will see later, this part will be important because it contains the search page result number.

The Craigslist website requires Javascript, so we'll need to use Selenium to load the page, and then use BeautifulSoup to extract the information we want.

2.1 Initialize a selenium driver and open Craigslist¶

As discussed in lecture, you can use Chrome, Firefox, or Edge as your selenium driver. In this part, you should do two things:

  1. Initialize the selenium driver
  2. Use the driver.get() function to open the following URL:

https://philadelphia.craigslist.org/search/apa?min_price=1&min_bedrooms=1&minSqft=1#search=1~gallery~0~0

This will give you the search results for 1-bedroom apartments in Philadelphia.

In [112]:
from selenium import webdriver
from selenium.webdriver.common.by import By
from selenium.webdriver.support.ui import WebDriverWait
from selenium.webdriver.support import expected_conditions as EC
from bs4 import BeautifulSoup
import pandas as pd
import requests
In [113]:
driver = webdriver.Chrome()
url = "https://philadelphia.craigslist.org/search/apa?min_price=1&min_bedrooms=1&minSqft=1#search=1gallery0~0"
driver.get(url)

2.2 Initialize your "soup"¶

Once selenium has the page open, we can get the page source from the driver and use BeautifulSoup to parse it. In this part, initialize a BeautifulSoup object with the driver's page source

In [114]:
aptsoup = BeautifulSoup(driver.page_source, "html.parser")

2.3 Parsing the HTML¶

Now that we have our "soup" object, we can use BeautifulSoup to extract out the elements we need:

  • Use the Web Inspector to identify the HTML element that holds the information on each apartment listing.
  • Use BeautifulSoup to extract these elements from the HTML.

At the end of this part, you should have a list of 120 elements, where each element is the listing for a specific apartment on the search page.

In [115]:
aptlistings = aptsoup.select(".cl-search-result")


len(aptlistings)
Out[115]:
120

2.4 Find the relevant pieces of information¶

We will now focus on the first element in the list of 120 apartments. Use the prettify() function to print out the HTML for this first element.

From this HTML, identify the HTML elements that hold:

  • The apartment price
  • The number of bedrooms
  • The square footage
  • The apartment title

For the first apartment, print out each of these pieces of information, using BeautifulSoup to select the proper elements.

In [116]:
aptlistings1 = aptlistings[0]
In [117]:
print(aptlistings1.prettify())

<li class="cl-search-result cl-search-view-mode-gallery" data-pid="7696705093" title="Fitness Center, Dishwasher, 1BD 1BA">
 <div class="gallery-card">
  <div class="cl-gallery">
   <div class="gallery-inner">
    <a class="main" href="https://philadelphia.craigslist.org/apa/d/philadelphia-fitness-center-dishwasher/7696705093.html">
     <div class="swipe" style="visibility: visible;">
      <div class="swipe-wrap" style="width: 5688px;">
       <div data-index="0" style="width: 316px; left: 0px; transition-duration: 0ms; transform: translateX(0px);">
        <span class="loading icom-">
        </span>
        <img alt="Fitness Center, Dishwasher, 1BD 1BA 1" src="https://images.craigslist.org/00i0i_VyzZhG8Sgq_09x06O_300x300.jpg"/>
       </div>
       <div data-index="1" style="width: 316px; left: -316px; transition-duration: 0ms; transform: translateX(316px);">
       </div>
       <div data-index="2" style="width: 316px; left: -632px; transition-duration: 0ms; transform: translateX(316px);">
       </div>
       <div data-index="3" style="width: 316px; left: -948px; transition-duration: 0ms; transform: translateX(316px);">
       </div>
       <div data-index="4" style="width: 316px; left: -1264px; transition-duration: 0ms; transform: translateX(316px);">
       </div>
       <div data-index="5" style="width: 316px; left: -1580px; transition-duration: 0ms; transform: translateX(316px);">
       </div>
       <div data-index="6" style="width: 316px; left: -1896px; transition-duration: 0ms; transform: translateX(316px);">
       </div>
       <div data-index="7" style="width: 316px; left: -2212px; transition-duration: 0ms; transform: translateX(316px);">
       </div>
       <div data-index="8" style="width: 316px; left: -2528px; transition-duration: 0ms; transform: translateX(-316px);">
       </div>
      </div>
     </div>
     <div class="slider-back-arrow icom-">
     </div>
     <div class="slider-forward-arrow icom-">
     </div>
    </a>
   </div>
   <div class="dots">
    <span class="dot selected">
     •
    </span>
    <span class="dot">
     •
    </span>
    <span class="dot">
     •
    </span>
    <span class="dot">
     •
    </span>
    <span class="dot">
     •
    </span>
    <span class="dot">
     •
    </span>
    <span class="dot">
     •
    </span>
    <span class="dot">
     •
    </span>
    <span class="dot">
     •
    </span>
   </div>
  </div>
  <a class="cl-app-anchor text-only posting-title" href="https://philadelphia.craigslist.org/apa/d/philadelphia-fitness-center-dishwasher/7696705093.html" tabindex="0">
   <span class="label">
    Fitness Center, Dishwasher, 1BD 1BA
   </span>
  </a>
  <div class="meta">
   12 mins ago
   <span class="separator">
    ·
   </span>
   <span class="housing-meta">
    <span class="post-bedrooms">
     1br
    </span>
    <span class="post-sqft">
     717ft
     <span class="exponent">
      2
     </span>
    </span>
   </span>
   <span class="separator">
    ·
   </span>
   1000 South Broad Street, Philadelphia, PA
  </div>
  <span class="priceinfo">
   $2,359
  </span>
  <button class="bd-button cl-favorite-button icon-only" tabindex="0" title="add to favorites list" type="button">
   <span class="icon icom-">
   </span>
   <span class="label">
   </span>
  </button>
  <button class="bd-button cl-banish-button icon-only" tabindex="0" title="hide posting" type="button">
   <span class="icon icom-">
   </span>
   <span class="label">
    hide
   </span>
  </button>
 </div>
</li>
In [118]:
price = aptlistings1.select_one(".priceinfo").text
price
Out[118]:
'$2,359'
In [119]:
numberofbeds=aptlistings1.select_one(".post-bedrooms").text
numberofbeds
Out[119]:
'1br'
In [120]:
sqft = aptlistings1.select_one(".post-sqft").text
sqft
Out[120]:
'717ft2'
In [121]:
apttitle =aptlistings1.select_one(".posting-title span.label").text
apttitle
Out[121]:
'Fitness Center, Dishwasher, 1BD 1BA'

2.5 Functions to format the results¶

In this section, you'll create functions that take in the raw string elements for price, size, and number of bedrooms and returns them formatted as numbers.

I've started the functions to format the values. You should finish theses functions in this section.

Hints

  • You can use string formatting functions like string.replace() and string.strip()
  • The int() and float() functions can convert strings to numbers
In [122]:
def format_bedrooms(bedrooms_string):
    nmbrbeds =''.join(n for n in bedrooms_string if n.isdigit())
    nmbrbeds = float(nmbrbeds)
    return nmbrbeds
In [123]:
def format_size(size_string):
    if size_string:
        size = ''.join(char for char in size_string if char.isdigit() or char in ('.', 'e', 'E', '+', '-'))
        return float(size)
In [124]:
def format_price(price_string):
    if price_string:
        price = ''.join(char for char in price_string if char.isdigit() or char in ('.', 'e', 'E', '+', '-'))
        return float(price)

2.6 Putting it all together¶

In this part, you'll complete the code block below using results from previous parts. The code will loop over 5 pages of search results and scrape data for 600 apartments.

We can get a specific page by changing the search=PAGE part of the URL hash. For example, to get page 2 instead of page 1, we will navigate to:

https://philadelphia.craigslist.org/search/apa?min_price=1&min_bedrooms=1&minSqft=1#search=2~gallery~0~0

In the code below, the outer for loop will loop over 5 pages of search results. The inner for loop will loop over the 120 apartments listed on each search page.

Fill in the missing pieces of the inner loop using the code from the previous section. We will be able to extract out the relevant pieces of info for each apartment.

After filling in the missing pieces and executing the code cell, you should have a Data Frame called results that holds the data for 600 apartment listings.

Notes¶

Be careful if you try to scrape more listings. Craigslist will temporarily ban your IP address (for a very short time) if you scrape too much at once. I've added a sleep() function to the for loop to wait 30 seconds between scraping requests.

If the for loop gets stuck at the "Processing page X..." step for more than a minute or so, your IP address is probably banned temporarily, and you'll have to wait a few minutes before trying again.

In [125]:
from time import sleep
In [126]:
from time import sleep
import pandas as pd
from bs4 import BeautifulSoup

results = []

# search in batches of 120 for 5 pages
max_pages = 5

# The base URL we will be using
base_url = "https://philadelphia.craigslist.org/search/apa?min_price=1&min_bedrooms=1&minSqft=1"

# loop over each page of search results
for page_num in range(1, max_pages + 1):
    print(f"Processing page {page_num}...")

    # Update the URL hash for this page number and make the combined URL
    url_hash = f"#search={page_num}~gallery~0~0"
    url = base_url + url_hash

    # Go to the driver and wait for 5 seconds
    driver.get(url)
    sleep(5)

    # YOUR CODE: get the list of all apartments
    # This is the same code from Part 1.2 and 1.3
    # It should be a list of 120 apartments
    
    soup = BeautifulSoup(driver.page_source, "html.parser")
    apts = soup.select(".cl-search-result")
    print("Number of apartments =", len(apts))

    # loop over each apartment in the list
    page_results = []
    for apt in apts:

        # YOUR CODE: the bedrooms string
        numberofbeds = apt.select_one(".post-bedrooms").text

        # YOUR CODE: the size string
        size = apt.select_one(".post-sqft").text

        # YOUR CODE: the title string
        title = apt.select_one(".posting-title span.label").text

        # YOUR CODE: the price string
        price = apt.select_one(".priceinfo").text

        # Format using functions from Part 1.5
        bedrooms = format_bedrooms(numberofbeds)
        size = format_size(size)
        price = format_price(price)

        # Save the result
        page_results.append([price, size, bedrooms, title])

    # Create a dataframe and save
    col_names = ["price", "size", "bedrooms", "title"]
    df = pd.DataFrame(page_results, columns=col_names)
    results.append(df)

    print("Sleeping for 10 seconds between calls")
    sleep(10)

# Finally, concatenate all the results
results = pd.concat(results, axis=0).reset_index(drop=True)

Processing page 1...
Number of apartments = 120
Sleeping for 10 seconds between calls
Processing page 2...
Number of apartments = 120
Sleeping for 10 seconds between calls
Processing page 3...
Number of apartments = 120
Sleeping for 10 seconds between calls
Processing page 4...
Number of apartments = 120
Sleeping for 10 seconds between calls
Processing page 5...
Number of apartments = 120
Sleeping for 10 seconds between calls

2.7 Plotting the distribution of prices¶

Use matplotlib's hist() function to make two histograms for:

  • Apartment prices
  • Apartment prices per square foot (price / size)

Make sure to add labels to the respective axes and a title describing the plot.

Side note: rental prices per sq. ft. from Craigslist¶

The histogram of price per sq ft should be centered around ~1.5. Here is a plot of how Philadelphia's rents compare to the other most populous cities:

No description has been provided for this image

Source

In [157]:
plt.hist(df['price'], bins=30, color='blue', edgecolor='black')
plt.title('Apartment prices')
plt.xlabel('Listing Price ($)')
plt.ylabel('Frequency')
plt.show()
No description has been provided for this image
In [158]:
plt.hist(df['price']/df['size'], bins=30, color='blue', edgecolor='black')
plt.title('Apartment prices per Square Foot')
plt.xlabel('Listing Price Per Square Foot($)')
plt.ylabel('Frequency')
plt.show()
No description has been provided for this image

2.8 Comparing prices for different sizes¶

Use altair to explore the relationship between price, size, and number of bedrooms. Make an interactive scatter plot of price (x-axis) vs. size (y-axis), with the points colored by the number of bedrooms.

Make sure the plot is interactive (zoom-able and pan-able) and add a tooltip with all of the columns in our scraped data frame.

With this sort of plot, you can quickly see the outlier apartments in terms of size and price.

In [129]:
 import altair as alt
In [130]:
import altair as alt
import pandas as pd

# Assuming your DataFrame is named results

# Scatter plot
scatter_plot = alt.Chart(results).mark_circle(size=60).encode(
    x=alt.X('price:Q', axis=alt.Axis(title='Price ($)')),
    y=alt.Y('size:Q', axis=alt.Axis(title='Size (sqft)')),
    color=alt.Color('bedrooms:N', scale=alt.Scale(scheme='category10')),
    tooltip=["title", "price:Q", "size:Q", "bedrooms:N"]
).properties(
    width=600,
    height=400
).interactive()

scatter_plot.display()
In [ ]: